**********************
*****mata program*****
**********************     
  
mata
version 11 
mata clear
mata set matastrict on

void  notch_mse_roundq (
            string scalar    incbin_name,
            string scalar    freq_name, 
			string scalar    yactual,
			real scalar      notch,       /* standardlized  */
            real scalar      low,     /* negative sign _number of bins below notch where excluded range for counterfactual starts */
            real scalar      high, 
            real scalar      degree,  /* degree of polynomial */
            real scalar      sf_s) 
{
real matrix         Poly,   get_beta_dummy, pick_right,  pick_bunch_and_right , dummies, x_nodum,  x, 
                    xxinv,  xxinvx,   x_plot, rounds , get_beta_polyround,get_beta_dummy_bunch,get_beta_dummy_hole,sround,
					get_beta_polyroundsround , get_beta_poly,get_beta_round,  get_beta_sround ,x_poly,x_poly_round ,
					get_beta_roundsround,x_roundsround 

        
real colvector     incbin, freq, y_actual, round80, round20, round50, iota,  cons, poly1, poly2, agg,  beta_hat,  beta_hat_dum, y_counter,
                   y_new,  y_hat,  y_hat_pred, error_adj, b_boot, m_boot, bn_boot,  bf_boot,  draw_j, resample,  freq_j,
                   y_counter_nosround, y_counter_noroundsround, y_actual_nosround, y_actual_noroundsround ,
				   beta_hat_j,  beta_hat_dum_j ,  beta_hat_polyround, beta_hat_dummy_bunch,beta_hat_dummy_hole, 
                   beta_hat_dummy_hole_j,beta_hat_dummy_bunch_j,Dummy,Z100, W50,
				   beta_hat_poly, beta_hat_round, beta_hat_sround,beta_hat_polyroundsround,
				   y_counter_j,               y_counter_nosround_j,  y_counter_noroundsround_j,
				   beta_hat_polyroundsround_j,beta_hat_polyround_j,   beta_hat_poly_j, y_roundsround,
				   beta_hat_roundsround
				          
real scalar       mn, rows, se,  last_poly_b, last_polyround_b, first_dummy_b, ndum,  xcols, minim,  maxim, maxim_rx, notch_rx, 
                  first_bunch_rx,  last_bunch_rx,  last_left_rx, first_right_rx,  people_nearbunch, bn, bf, b, 
			       n_it , bn_j, bf_j, b_j, mp_j, m_j, n_it_j, i,j,ii, jj, b_se, m_se, nround, ndumbunch,ndumhole,first_hole_rx, 
				   last_hole_rx,mp,af, m, af_boot,af_j,nsround, diff, tolerance, mse 				   
				  
string scalar    b_s		
			  
/* GETTING MATRICES FOR ESTIMATION */	
  incbin			=st_data(.,incbin_name) 
  freq			    =st_data(.,freq_name)
  y_actual          =st_data(.,yactual)
 
  round50           =st_data(.,"round50")
  Z100              =st_data(.,"Z100")
  rounds            =round50, Z100

/* CONSTRUCTING CONSTANT, POLYNOMIAL TERM, AND DUMMIES TERM */
  rows            =rows(incbin)
  iota            =J(rows,1,1)
  cons            =iota			
			
/* NORMALIZE PLOY1 - DEMEAN AND DIVIDE BY SE   */
  mn            =(iota'*incbin)/(iota'*iota)    
  poly1         =incbin:-mn                 
  se            =(sqrt(cross(poly1,poly1)/rows)) 
  poly1         =poly1/se                    
  Poly          =poly1
  for (i=2; i<=degree; i++) Poly = Poly,poly1:^i	
  poly2         =poly1:^2	
  sround=Z100:*poly1, Z100:*poly2, round50:*poly1, round50:*poly2

/* ROW NUMBERS IN BETA */
  nround          =2
  nsround         =4
  ndumhole        =-low	         //# of hole dummies
  ndumbunch       =high+1
  last_polyround_b=1+degree+nround+nsround
  first_dummy_b   =1+degree+nround+nsround+1  

/*defining notch and bunching window */	
  ndum  = high-low+1             //# of dummies
  xcols = 1+degree+nround+nsround+ndum 

/* ROW NUMBERS FOR BUNCHING RANGE */
  minim = colmin(incbin) 		 /* lowest incbin  */
  maxim = colmax(incbin)         /* highest incbin  */   
			
/* ROW NUMNBERS FOR EXCLUDED RANGE */	
  maxim_rx        = (maxim-minim)/sf_s+1     
  notch_rx        =  notch-minim /sf_s+1     /*bin # of notch*/       
	
  first_hole_rx  = notch_rx+low	             /*bin # of first_bunch, usually 1 if minim=notch+low*/
  last_hole_rx   = notch_rx-1
  first_bunch_rx = notch_rx
  last_bunch_rx  = notch_rx+high		   
  last_left_rx   = first_hole_rx-1
  first_right_rx = last_bunch_rx+1

/* USEFUL SCALARS, VECTORS AND MATRICES */
  get_beta_polyroundsround      = diag(J((1+degree+nround+nsround ),1,1))  , J((1+degree+nround+nsround ),ndum,0)
  get_beta_dummy_hole           = J(ndumhole,(1+degree+nround+nsround),0) , diag(J(ndumhole,1,1))  , J(ndumhole,ndumbunch,0)
  get_beta_dummy_bunch          = J(ndumbunch,(1+degree+nround+nsround),0) ,J(ndumbunch,ndumhole,0), diag(J(ndumbunch,1,1))					 
  get_beta_roundsround          = J((nround+nsround),(1+degree),0)         , diag(J((nround+nsround),1,1)),  J((nround+nsround),ndum,0)  

  /* CREATING DUMMIES FOR EXCLUDED RANGE IN COUNTERFACTUAL */
  dummies = J(rows, ndum, .) 
  ii=0
  for(jj=first_hole_rx; jj<=last_bunch_rx; jj++) {
        ii = ii+1
       dummies[.,ii] = e(jj,rows)'
   }
		      
/* CREATING THE MATRICES NEEDED FOR ESTIMATION OF BETA_HAT */
  x_nodum = cons, Poly, rounds, sround
  x = x_nodum, dummies 
  xxinv = invsym(cross(x,x))
  xxinvx = xxinv*x'
  
  x_poly=cons, Poly
  x_poly_round=cons, Poly, rounds
  x_roundsround=rounds, sround
 /* CONSTRUCTING BETAHAT AND Y_COUNTER*/
  beta_hat                 =xxinvx*freq 
  beta_hat_polyroundsround =get_beta_polyroundsround *beta_hat	
  beta_hat_dummy_hole      =get_beta_dummy_hole      *beta_hat
  beta_hat_dummy_bunch     =get_beta_dummy_bunch     *beta_hat
  
  beta_hat_roundsround     =get_beta_roundsround     *beta_hat

	
 /* CONSTRUCTING FITTED DATA */ 
  y_counter =x_nodum*beta_hat_polyroundsround		
  y_hat     =x*beta_hat

  y_roundsround =x_roundsround*beta_hat_roundsround 

 /* CALCULATING EXCESS MASS AND ELASTICITY/EVASION RESPONSE */
  mp = colsum(beta_hat_dummy_hole)
  bn = colsum(beta_hat_dummy_bunch)
  diff=abs(bn+mp)
  if (bn>mp)  tolerance=bn*0.1
  if (bn<mp)  tolerance=mp*0.1
  
if (diff<=tolerance) {	
	 mse=((y_actual-y_hat/4)'*(y_actual-y_hat/4))/rows	  
	 }
     else {
	    mse=. 
		}
	
	st_global("mse", strofreal(mse,"%9.2f"))			
   }
mata mosave  notch_mse_roundq(),replace
end
   
   
